#include <iostream>
#include <math.h>
#include <string.h>
#include "cmes_springxcloud_business_sysIdentify_util_WindStream.h"
using namespace std;


int g_arrowden =  20 ;

typedef struct S_Strmpoint
{
  int lineState; //线条起点还是续点--起点为0 还是续点为非0，一般设为1
  float x,y;
  int color;
  struct S_Strmpoint * next;
}S_STRMPOINT;

/*记录点的结构 可以记录普通点，也可记录箭头--记录两个点和一个颜色*/
/*线条记录一个点，一个颜色，一个起点还是续点的标记*/
typedef struct S_Arrow
{
  float x,y;
  float x2,y2;
  int color;
  struct S_Arrow * next;
}S_ARROW;

static int intflg;                        /* Batch flag */
static int color;                         /* Current color      */

//static float a150=150.0*3.1416/180.0;
float StartLon,EndLon,StartLat,EndLat,LonInterval,LatInterval;

S_ARROW *arrow=NULL;
S_STRMPOINT *streamline=NULL;
S_ARROW *arrowT, *arrowR;
S_STRMPOINT *streamlineT,*streamlineR;

// 函数声明 ----------------------------------------------------------------------------------
void strmShadFree();
float *write2Array(float *u, float *v, int isiz, int jsiz);
void writeLineOut(const char *filename, float *array);
float *readWindData(const char *file, int den);


// 函数声明 end ------------------------------------------------------------------------------

/*
//20090330 流线分析--增加参数，适合各类数据分析，主要是增加对等距格点的支持
extern "C" float *  streamAnalyse_f (float *u, float *v, float *c, int is, int js,
   float uund, float vund, float cund, int flag, float *shdlvs,
   int *shdcls, int shdcnt, int den,float SLon,float SLat,float XInterval ,float YInterval)


   //使用WINAPI参数，可以编译通用的DLL,可以由其它语言，如VB,FORTRAN等调用,但要包含windows.h文件
extern "C" float *  streamAnalyse (float *u, float *v, float *c, int is, int js,
      float uund, float vund, float cund, int flag, float *shdlvs,
     int *shdcls, int shdcnt, int den);

*/

// main函数用于测试。
// int main()
// {
//     cout << "Hello world!" << endl;
//     float * retArray;
//     int den = 10 ;//流线密度
//     retArray=readWindData("E:\\coding\\SMM-Projects\\stream\\22062600.012",den );
//     writeLineOut("22062600.012.output24.000.txt",retArray);
//     free(retArray);
//     printf("press any key to quit.\n") ;
//     getchar() ;
//     return 0;
// }

JNIEXPORT void JNICALL Java_cmes_springxcloud_business_sysIdentify_util_WindStream_test
  (JNIEnv *, jobject )
{
	cout<<"test test test"<<endl;
}

/**
 * 输入micaps数据生成风流线结果文件。
 * @params infilenames 输入文件完整路径。
 * @params grid_density 流线密度，取值1-10，或者大于10.
 * @params outfilename 输出结果完整文件路径。
 * @return 成功返回0，失败返回非零数字。
 *
 */
JNIEXPORT jint JNICALL Java_cmes_springxcloud_business_sysIdentify_util_WindStream_run
  (JNIEnv* env, jobject jobj, jstring jinfilename, jint jden,jint jarrowden, jstring joutfilename)
{
	cout<<"Java_cmes_springxcloud_business_sysIdentify_util_WindStream_run start..."<<endl;
	int den = jden ;
	const char* infilename = (env)->GetStringUTFChars(  jinfilename, NULL ) ;
	const char* outfilename = (env)->GetStringUTFChars(  joutfilename, NULL ) ;
	cout<<"infilename: "<<infilename<<endl;
	cout<<"outfilename: "<<outfilename<<endl;
	cout<<"den: "<<den<<endl;
	cout<<"arrow den: "<<jarrowden<<endl;
	
	g_arrowden = jarrowden ;

	float* retArray=readWindData(infilename ,den );
	writeLineOut(outfilename,retArray);
	free(retArray);

	cout<<"Java_cmes_springxcloud_business_sysIdentify_util_WindStream_run end."<<endl;
	return 0 ;
}


void gxplot(float xx, float yy, int cmd)
{
	// printf ("%f %f %d %d\n",xx,yy,cmd,color);
	static int _s_debug_gxplot_streamline = 0 ;
	++_s_debug_gxplot_streamline ;
	//if(_s_debug_gxplot_streamline%10000==0)printf("_s_debug_gxplot_streamline %d\n" , _s_debug_gxplot_streamline);

	streamlineT = (S_STRMPOINT *)malloc(sizeof(S_STRMPOINT));
	streamlineT->color = color;
	streamlineT->x = xx;
	streamlineT->y = yy;
	streamlineT->next = NULL;
	if (cmd == 3)
		streamlineT->lineState = 0;
	else if (cmd == 2)
		streamlineT->lineState = 1;
	else
		streamlineT->lineState = -1;

	if (streamline == NULL)
	{
		streamline = streamlineT;
		streamlineR = streamline;
	}
	else
	{
		streamlineR->next = streamlineT;
		streamlineR = streamlineT;
	}
}
void gxcolr(int clr)
{ /* Set color     */
	if (clr < 0)
		clr = 0;
	if (clr > 99)
		clr = 99;
	color = clr;
}


void strmar(float xx1, float yy1, float xx2, float yy2)
{
    static int _s_debug_strmar = 0 ;
    ++_s_debug_strmar;
    //if( _s_debug_strmar%10000==0 ) {
    //    printf("_s_debug_strmar %d\n" , _s_debug_strmar) ;
    //}


	arrowT = (S_ARROW *)malloc(sizeof(S_ARROW));
	arrowT->color = color;
	arrowT->x = xx1;
	arrowT->y = yy1;
	arrowT->x2 = xx2;
	arrowT->y2 = yy2;
	arrowT->next = NULL;

	if (arrow == NULL)
	{
		arrow = arrowT;
		arrowR = arrow;
		/*printf("��һ��arrow %f %f\n",xx1,yy1);*/
	}
	else
	{
		/*printf("      arrow %f %f\n",xx1,yy1);*/
		arrowR->next = arrowT;
		arrowR = arrowT;
	}
}
int gxshdc(float *shdlvs, int *shdcls, int shdcnt, float val)
{
	int i;

	if (shdcnt == 0)
		return (1);
	if (shdcnt == 1)
		return (shdcls[0]);
	if (val < shdlvs[1])
		return (shdcls[0]);
	for (i = 1; i < shdcnt - 1; i++)
	{
		if (val >= shdlvs[i] && val < shdlvs[i + 1])
			return (shdcls[i]);
	}
	return (shdcls[shdcnt - 1]);
}

void gxconv(float s, float t, float *x, float *y) //, int level)
{
	*x = (StartLon + (s - 1) * LonInterval);
	*y = (StartLat + (t - 1) * LatInterval);
}

float *streamAnalyse(float *u, float *v, float *c, // u,v分量   c风速
                     int is, int js, // 网格尺寸
								float uund, float vund, float cund, // u，v，c无效值
								int flag, //是否使用颜色
								float *shdlvs, //颜色映射的值
								int *shdcls,   //对应颜色
								int shdcnt,    //颜色数量
								int den        //流线密度
								)
{
	float *up, *vp;
	float x, y, xx, yy, uv1, uv2, uv, vv1, vv2, vv, auv, avv, xsav, ysav, xold, yold;
	int i, ii, jj, ii1, ij1, i2, j2, ipt, acnt, icol, scol, dis; // kk,ipen,j
	int *it, siz, iacc, iisav, iscl, imn, imx, jmn, jmx, iz, jz, iss, jss, bflg;
	float fact, rscl, xxsv, yysv, *cp, cv1, cv2, cv;

	printf("streamAnalyse is %d js %d shdcnt %d\n" , is,js,den) ;

	float *retArray;

	scol = -9;
	icol = 1;

	/* Figure out the interval for the flag grid */

	i = is;
	if (js > i)
		i = js;
	iscl = 200 / i;
	iscl = iscl + den - 5;
	if (iscl < 1)
		iscl = 1;
	if (iscl > 10)
		iscl = 10;
	fact = (float)0.5 / ((float)iscl);
	rscl = (float)iscl;

	/* Allocate memory for the flag grid */

	iss = is * iscl;
	jss = js * iscl;
	siz = iss * jss;
	printf("iss %d jss %d\n",iss,jss) ;
	it = (int *)malloc(sizeof(int) * siz);
	if (it == NULL)
	{
		printf("Cannot allocate memory for streamline function\n");
		return NULL;
	}
	for (i = 0; i < siz; i++)
		*(it + i) = 0;

	/* Loop through flag grid to look for start of streamlines.
	   To start requires no streams drawn within surrounding
	   flag boxes.  */

	i2 = 0;
	j2 = 0;
	for (i = 0; i < siz; i++)
	{
		dis = 2;
		if (den < 5)
			dis = 3;
		if (den > 5)
			dis = 1;
		imn = i2 - dis;
		imx = i2 + dis + 1;
		jmn = j2 - dis;
		jmx = j2 + dis + 1;
		if (imn < 0)
			imn = 0;
		if (imx > iss)
			imx = iss;
		if (jmn < 0)
			jmn = 0;
		if (jmx > jss)
			jmx = jss;
		iacc = 0;
		for (jz = jmn; jz < jmx; jz++)
		{
			ipt = jz * iss + imn;
			for (iz = imn; iz < imx; iz++)
			{
				iacc = iacc + *(it + ipt);
				ipt++;
			}
		}
		if (iacc == 0)
		{
			x = ((float)i2) / rscl;
			y = ((float)j2) / rscl;
			xsav = x;
			ysav = y;
			gxconv(x + 1.0, y + 1.0, &xx, &yy); //,3);
			gxplot(xx, yy, 3);
			xxsv = xx;
			yysv = yy;
			iisav = -999;
			iacc = 0;
			acnt = 0;
			bflg = 0;
			while (x >= 0.0 && x < (float)(is - 1) && y >= 0.0 && y < (float)(js - 1))
			{
				ii = (int)x;
				jj = (int)y;
				xx = x - (float)ii;
				yy = y - (float)jj;
				up = u + jj * is + ii;
				vp = v + jj * is + ii;
				if (*up == uund || *(up + 1) == uund ||
					*(up + is) == uund || *(up + is + 1) == uund)
					break;
				if (*vp == vund || *(vp + 1) == vund ||
					*(vp + is) == vund || *(vp + is + 1) == vund)
					break;
				if (flag)
				{
					cp = c + jj * is + ii;
					if (*cp == cund || *(cp + 1) == cund ||
						*(cp + is) == cund || *(cp + is + 1) == cund)
						icol = 15;
					else
					{
						cv1 = *cp + (*(cp + 1) - *cp) * xx;
						cv2 = *(cp + is) + (*(cp + is + 1) - *(cp + is)) * xx;
						cv = cv1 + (cv2 - cv1) * yy;
						icol = gxshdc(shdlvs, shdcls, shdcnt, cv);
					}
					if (icol != scol && icol > -1)
						gxcolr(icol);
					scol = icol;
				}
				uv1 = *up + (*(up + 1) - *up) * xx;
				uv2 = *(up + is) + (*(up + is + 1) - *(up + is)) * xx;
				uv = uv1 + (uv2 - uv1) * yy;
				vv1 = *vp + (*(vp + 1) - *vp) * xx;
				vv2 = *(vp + is) + (*(vp + is + 1) - *(vp + is)) * xx;
				vv = vv1 + (vv2 - vv1) * yy;
				auv = fabs(uv);
				avv = fabs(vv);
				if (auv < 0.1 && avv < 0.1)
					break;
				if (auv > avv)
				{
					vv = vv * fact / auv;
					uv = uv * fact / auv;
				}
				else
				{
					uv = uv * fact / avv;
					vv = vv * fact / avv;
				}
				x = x + uv;
				y = y + vv;
				ii1 = (int)(x * rscl);
				ij1 = (int)(y * rscl);
				ii1 = ij1 * iss + ii1;
				if (ii1 < 0 || ii1 >= siz)
					break;
				if (*(it + ii1) == 1)
					break;
				if (ii1 != iisav && iisav > -1)
					*(it + iisav) = 1;
				if (ii1 == iisav)
					iacc++;
				else
					iacc = 0;
				if (iacc > 10)
					break;
				iisav = ii1;
				gxconv(x + 1.0, y + 1.0, &xx, &yy); //,3);
				if (icol > -1)
				{
					if (bflg)
					{
						gxplot(xold, yold, 3);
						bflg = 0;
					}
					gxplot(xx, yy, 2);
				}
				else
					bflg = 1;
				xold = xx;
				yold = yy;
				acnt++;
				if (acnt > 20 )
				{
					if (icol > -1)
						strmar(xxsv, yysv, xx, yy);
					acnt = 0;
				}
				xxsv = xx;
				yysv = yy;
			}
			bflg = 0;
			x = xsav;
			y = ysav;
			gxconv(x + 1.0, y + 1.0, &xx, &yy); //,3);
			gxplot(xx, yy, 3);
			xxsv = xx;
			yysv = yy;
			iisav = -999;
			iacc = 0;
			acnt = 19;
			while (x >= 0.0 && x < (float)(is - 1) && y >= 0.0 && y < (float)(js - 1))
			{
				ii = (int)x;
				jj = (int)y;
				xx = x - (float)ii;
				yy = y - (float)jj;
				up = u + jj * is + ii;
				vp = v + jj * is + ii;
				if (*up == uund || *(up + 1) == uund ||
					*(up + is) == uund || *(up + is + 1) == uund)
					break;
				if (*vp == vund || *(vp + 1) == vund ||
					*(vp + is) == vund || *(vp + is + 1) == vund)
					break;
				if (flag)
				{
					cp = c + jj * is + ii;
					if (*cp == cund || *(cp + 1) == cund ||
						*(cp + is) == cund || *(cp + is + 1) == cund)
						icol = 15;
					else
					{
						cv1 = *cp + (*(cp + 1) - *cp) * xx;
						cv2 = *(cp + is) + (*(cp + is + 1) - *(cp + is)) * xx;
						cv = cv1 + (cv2 - cv1) * yy;
						icol = gxshdc(shdlvs, shdcls, shdcnt, cv);
					}
					if (icol != scol && icol > -1)
						gxcolr(icol);
					scol = icol;
				}
				uv1 = *up + (*(up + 1) - *up) * xx;
				uv2 = *(up + is) + (*(up + is + 1) - *(up + is)) * xx;
				uv = uv1 + (uv2 - uv1) * yy;
				vv1 = *vp + (*(vp + 1) - *vp) * xx;
				vv2 = *(vp + is) + (*(vp + is + 1) - *(vp + is)) * xx;
				vv = vv1 + (vv2 - vv1) * yy;
				auv = fabs(uv);
				avv = fabs(vv);
				if (auv < 0.1 && avv < 0.1)
					break;
				if (auv > avv)
				{
					vv = vv * fact / auv;
					uv = uv * fact / auv;
				}
				else
				{
					uv = uv * fact / avv;
					vv = vv * fact / avv;
				}
				x = x - uv;
				y = y - vv;
				ii1 = (int)(x * rscl);
				ij1 = (int)(y * rscl);
				ii1 = ij1 * iss + ii1;
				if (ii1 < 0 || ii1 >= siz)
					break;
				if (*(it + ii1) == 1)
					break;
				if (ii1 != iisav && iisav > -1)
					*(it + iisav) = 1;
				if (ii1 == iisav)
					iacc++;
				else
					iacc = 0;
				if (iacc > 10)
					break;
				iisav = ii1;
				gxconv(x + 1.0, y + 1.0, &xx, &yy); //,3);
				if (icol > -1)
				{
					if (bflg)
					{
						gxplot(xold, yold, 3);
						bflg = 0;
					}
					gxplot(xx, yy, 2);
				}
				else
					bflg = 1;
				xold = xx;
				yold = yy;
				acnt++;
				if (acnt > 20)
				{
					if (icol > -1)
						strmar(xx, yy, xxsv, yysv);
					acnt = 0;
				}
				xxsv = xx;
				yysv = yy;
			}
		}
		i2++;
		if (i2 == iss)
		{
			i2 = 0;
			j2++;
		}
	}
	free(it);
	retArray = write2Array(u, v, is, js);
	strmShadFree();
	return retArray;
}

float *readWindData(const char *file, int den)
{
	FILE *fuv;
	char title[10];
	char Chart_Name[100];
	// char filename[200];
	int idata_type;
	int iYear, iMonth, iDay, iTime, iValid, iLevel;
	float iIntervalX, iIntervalY;
	int iSize, jSize;
	// int inFileHandle;
	int i, j;
	float *u;
	float *v;
	float *c;
	float tmplat;
	// int tmpi;
	float uund, vund, cund;
	// int den;
	int flag;
	float *shdlvs;
	int *shdcls;
	int shdcnt;
	// int DefColor;

	int ret;
	float *retArray;
	// Note: file would be modified by system (add one character "8")
	//       when apply new memory u(isize*jsize) on Redhat 9.
	// strcpy(filename,file);
	retArray = (float *)malloc(sizeof(float));

	if ((fuv = fopen(file, "r")) == NULL)
	{
		fprintf(stderr, "Cannot open wind data file.\n");
		retArray[0] = -1;
		return retArray;
	}
	ret = fscanf(fuv, "%s%d%s\n", title, &idata_type, Chart_Name);
	if (ret != 3)
	{
		retArray[0] = -2;
		return retArray;
	}
	ret = strcmp(title, "diamond");
	if (ret != 0)
	{
		retArray[0] = -3;
		return retArray;
	}

	ret = fscanf(fuv, "%d %d %d %d %d %d %f %f %f %f %f %f %d %d\n", &iYear, &iMonth, &iDay,
				 &iTime, &iValid, &iLevel, &iIntervalX, &iIntervalY, &StartLon, &EndLon, &StartLat,
				 &EndLat, &iSize, &jSize);
	if (ret != 14)
	{
		retArray[0] = -4;
		return retArray;
	}

	u = (float *)malloc(iSize * jSize * sizeof(float));
	v = (float *)malloc(iSize * jSize * sizeof(float));
	c = (float *)malloc(iSize * jSize * sizeof(float));

	if (iIntervalY > 0)
	{
		for (j = 0; j < jSize * iSize; j++)
		{
			ret = fscanf(fuv, "%f", &u[j]);
			if (ret != 1)
			{
				retArray[0] = -5;
				return retArray;
			}
		}
		for (j = 0; j < jSize * iSize; j++)
		{
			ret = fscanf(fuv, "%f", &v[j]);
			if (ret != 1)
			{
				retArray[0] = -6;
				return retArray;
			}
		}
	}
	else
	{
		iIntervalY = -iIntervalY;
		tmplat = StartLat;
		StartLat = EndLat;
		EndLat = tmplat;

		for (j = jSize - 1; j >= 0; j--)
		{
			for (i = 0; i < iSize; i++)
			{
				// tmpi=i+j*iSize;
				ret = fscanf(fuv, "%f", &u[i + j * iSize]);

				if (ret != 1)
				{
					retArray[0] = -7;
					return retArray;
				}
			}
		}

		for (j = jSize - 1; j >= 0; j--)
		{
			for (i = 0; i < iSize; i++)
			{

				ret = fscanf(fuv, "%f", &v[i + j * iSize]);

				if (ret != 1)
				{
					retArray[0] = -8;
					return retArray;
				}
			}
		}
	}
	fclose(fuv);

	free(retArray);

	LonInterval = iIntervalX;
	LatInterval = iIntervalY;
	for (j = jSize - 1; j >= 0; j--)
	{
		for (i = 0; i < iSize; i++)
		{
			c[i + j * iSize] = sqrt(u[i + j * iSize] * u[i + j * iSize] + v[i + j * iSize] * v[i + j * iSize]);
		}
	}

	shdcnt = 40;
	flag = 1;
	uund = 9999.0;
	vund = 9999.0;
	cund = 9999.0;
	// den=5;

	shdcls = (int *)malloc(shdcnt * sizeof(int));
	shdlvs = (float *)malloc(shdcnt * sizeof(float));
	for (i = 0; i < 40; i++)
	{
		shdcls[i] = i;
		shdlvs[i] = (float)i * 2;
	}

	retArray = streamAnalyse(u, v, c, iSize, jSize, uund, vund, cund, flag,
							 shdlvs, shdcls, shdcnt, den);
	free(u);
	free(v);
	free(c);
	free(shdlvs);
	free(shdcls);
	return retArray;
}

void strmShadFree()
{
	if (streamline != NULL)
	{
		while (streamline != NULL)
		{
			streamlineT = streamline;
			streamline = streamline->next;
			free(streamlineT);
		} // stream line
	}
	while (arrow != NULL)
	{
		arrowT = arrow;
		arrow = arrow->next;
		free(arrowT);
	}
}

float *write2Array(float *u, float *v, int isiz, int jsiz)
{
	float *retArray;
	int i, j;
	long iorder;
	int arrowNum;
	int arraySize;
	arraySize = 0;
	/*printf("��������� \n");*/
	if (streamline == NULL)
	{
		retArray = NULL;
	}
	else
	{
		// ��������ͼ�ͷ
		/*printf("���������...\n");*/
		streamlineT = streamline;
		while (streamlineT != NULL)
		{
			// First point, just want to minimize if-then-else
			arraySize++;
			streamlineT = streamlineT->next;
		} // stream line
		printf("streamline arraySize %d\n" , arraySize) ;
		/*printf("    ��ռ�������С %d\n",arraySize*4);*/
		if (arrow == NULL)
		{
			printf("    û�м�ͷ \n");
		}
		arrowT = arrow;
		arrowNum = 0;
		while (arrowT != NULL)
		{
			arrowNum++;
			arrowT = arrowT->next;
		}
		printf("arrowNum %d\n",arrowNum);
		/*printf("    ��ͷ���� %d ռ���ڴ�%d\n",arrowNum,arrowNum*6);*/
		int lineSize = arraySize * 4 + arrowNum * 6 + 2;
		arraySize = arraySize * 4 + arrowNum * 6 + 7 + isiz * jsiz * 2 + 4;
		/*printf("    �����С %d\n",arraySize);*/
		retArray = (float *)malloc(arraySize * sizeof(float));

		retArray[0] = arraySize;
		retArray[1] = lineSize;

		arraySize = 0;
		streamlineT = streamline;
		while (streamlineT != NULL)
		{
			// First point, just want to minimize if-then-else
			retArray[2 + arraySize * 4] = streamlineT->lineState;
			retArray[2 + arraySize * 4 + 1] = streamlineT->x;
			retArray[2 + arraySize * 4 + 2] = streamlineT->y;
			retArray[2 + arraySize * 4 + 3] = streamlineT->color;
			// printf(" %d %d %f %f\n",arraySize*4,streamlineT->lineState,streamlineT->x,streamlineT->y);
			arraySize++;
			streamlineT = streamlineT->next;
		} // stream line
        printf("arraySize2 %d\n",arraySize);
		// draw arrow
		arrowT = arrow;
		arrowNum = 0;
		while (arrowT != NULL)
		{
			// u = arrow->x1;
			// v = arrow->y1;
			// u2= arrow->x2;
			// v2= arrow->y2;
			// dir = atan2(y-y2,x-x2);

			retArray[2 + arraySize * 4 + arrowNum * 6] = 2;
			retArray[2 + arraySize * 4 + arrowNum * 6 + 1] = arrowT->x;
			retArray[2 + arraySize * 4 + arrowNum * 6 + 2] = arrowT->y;
			retArray[2 + arraySize * 4 + arrowNum * 6 + 3] = arrowT->x2;
			retArray[2 + arraySize * 4 + arrowNum * 6 + 4] = arrowT->y2;
			retArray[2 + arraySize * 4 + arrowNum * 6 + 5] = arrowT->color; // �Ƕ�û��ͶӰ����ͼ��
			// printf(" %d %d %f %f\n",arraySize*4+arrowNum*6,arrowT->color,arrowT->x,arrowT->y);
			arrowNum++;
			arrowT = arrowT->next;
		}
		printf("arrowNum2 %d\n",arrowNum);
		// retArray[arraySize*4+arrowNum*6+5]=-1;
	}
	/*printf("��������\n");*/
	iorder = 2 + arraySize * 4 + arrowNum * 6;
	iorder++;
	retArray[iorder] = -1;
	iorder++;
	retArray[iorder] = isiz;
	iorder++;
	retArray[iorder] = jsiz;
	iorder++;
	retArray[iorder] = StartLon;
	iorder++;
	retArray[iorder] = StartLat;
	iorder++;
	retArray[iorder] = EndLon;
	iorder++;
	retArray[iorder] = EndLat;
	iorder++;
	retArray[iorder] = LonInterval;
	iorder++;
	retArray[iorder] = LatInterval;
	for (j = 0; j < jsiz; j++)
	{
		for (i = 0; i < isiz; i++)
		{
			iorder++;
			retArray[iorder] = u[i + j * isiz];
			iorder++;
			retArray[iorder] = v[i + j * isiz];
		}
	}
	return retArray;
}
void writeLineOut(const char *filename, float *array)
{
	FILE *fout;
	int i, k;
	float a, b, c, d, e, f, g, h;
	int arrowCursor = 0 ;
	/*printf("����������ļ�\n");*/
	if (array == NULL)
	{
		// printf("������\n");
		printf("error: array is null.\n");//by wf
	}
	if ((fout = fopen(filename, "w")) == NULL)
	{
		fprintf(stderr, "Cannot open output data file\n");
		return;
	}
	i = 0;
	a = array[i];
	if ((int)a == 0)
	{
		fclose(fout);
		return;
	}
	i++;
	i++;
	while (i >= 0)
	{
		a = array[i];
		i++;
		if ((int)(a) == 1 || (int)(a) == 0)
		{
			c = array[i];
			i++;
			d = array[i];
			i++;
			e = array[i];
			i++;
			// a 线端点标志（0-起点，1-中间或者终点）；c-点经度坐标；d-点纬度坐标；e-线颜色标志。
			fprintf(fout, "%f %f %f %f\n", a, c, d, e);
		}
		else if ((int)(a) == 2)
		{
			c = array[i];
			i++;
			d = array[i];
			i++;
			e = array[i];
			i++;
			f = array[i];
			i++;
			g = array[i];
			i++;
			++ arrowCursor ;
			if( arrowCursor%g_arrowden==0 ){
				fprintf(fout, "%f %f %f %f %f %f\n", a, c, d, e, f, g);
			}
			
		}
		else
		{
			// print u and v value out
			a = array[i];
			i++;
			fprintf(fout,"线条结束标志 %f\n",a);
			// wf fprintf(fout, "����������־ %f\n", a);
			a = array[i];
			i++;
			b = array[i];
			i++;
			c = array[i];
			i++;
			d = array[i];
			i++;
			e = array[i];
			i++;
			f = array[i];
			i++;
			g = array[i];
			i++;
			h = array[i];
			i++;
			fprintf(fout, "%f %f %f %f %f %f %f %f\n", a, b, c, d, e, f, g, h);
			for (k = 0; k < a * b; k++)
			{
				c = array[i];
				i++;
				d = array[i];
				i++;
				fprintf(fout, "%f %f\n", c, d);
			}
			//}
			fclose(fout);
			return;
		}
	}

	return;
}

//20090330 流线分析--增加参数，适合各类数据分析，主要是增加对等距格点的支持
float *  streamAnalyse_f (float *u, float *v, float *c, int is, int js,
   float uund, float vund, float cund, int flag, float *shdlvs,
   int *shdcls, int shdcnt, int den,float SLon,float SLat,float XInterval ,float YInterval)
{


	streamline=NULL; //20090403 lya
	arrow=NULL;

  float *up, *vp;
  float x,y,xx,yy,uv1,uv2,uv,vv1,vv2,vv,auv,avv,xsav,ysav,xold,yold;
  int i,ii,jj,ii1,ij1,i2,j2,ipt,acnt,icol,scol,dis;//kk,ipen,j
  int *it,siz,iacc,iisav,iscl,imn,imx,jmn,jmx,iz,jz,iss,jss,bflg;
  float fact,rscl,xxsv,yysv,*cp,cv1,cv2,cv;

  float * retArray;

  scol = -9;
  icol = 1;

  /* Figure out the interval for the flag grid */

  i = is;
  if (js>i) i = js;
  iscl = 200/i;
  iscl = iscl + den - 5;
  if (iscl<1) iscl=1;
  if (iscl>10) iscl=10;
  fact = (float)0.5/((float)iscl);
  rscl = (float)iscl;

  /* Allocate memory for the flag grid */

  iss = is*iscl; jss = js*iscl;
  siz = iss*jss;
  it = (int *)malloc(sizeof(int) * siz);
  if (it==NULL) {
    //printf ("Cannot allocate memory for streamline function\n");
    return NULL;
  }
  for (i=0; i<siz; i++) *(it+i) = 0;

  /* Loop through flag grid to look for start of streamlines.
     To start requires no streams drawn within surrounding
     flag boxes.  */

  i2 = 0;
  j2 = 0;
  for (i=0; i<siz; i++) {
    dis = 2;
    if (den<5) dis = 3;
    if (den>5) dis = 1;
    imn = i2-dis; imx = i2+dis+1;
    jmn = j2-dis; jmx = j2+dis+1;
    if (imn<0) imn = 0;
    if (imx>iss) imx = iss;
    if (jmn<0) jmn = 0;
    if (jmx>jss) jmx = jss;
    iacc = 0;
    for (jz=jmn; jz<jmx; jz++) {
      ipt = jz*iss+imn;
      for (iz=imn; iz<imx; iz++) {
        iacc = iacc + *(it+ipt);
        ipt++;
      }
    }
    if (iacc==0){
      x = ((float)i2)/rscl;
      y = ((float)j2)/rscl;
      xsav = x;
      ysav = y;
      // wf -----------------------   gxconv_lonlat (x,y,&xx,&yy,SLon, SLat, XInterval , YInterval);//,3);
      gxplot (xx,yy,3);
      xxsv = xx; yysv = yy;
      iisav = -999;
      iacc = 0;
      acnt = 0;
      bflg = 0;
      while (x>=0.0 && x<(float)(is-1) && y>=0.0 && y<(float)(js-1)) {
        ii = (int)x;
        jj = (int)y;
        xx = x - (float)ii;
        yy = y - (float)jj;
        up = u + jj*is+ii;
        vp = v + jj*is+ii;
        if (*up==uund || *(up+1)==uund ||
            *(up+is)==uund || *(up+is+1)==uund) break;
        if (*vp==vund || *(vp+1)==vund ||
            *(vp+is)==vund || *(vp+is+1)==vund) break;
        if (flag) {
          cp = c + jj*is+ii;
          if (*cp==cund || *(cp+1)==cund ||
              *(cp+is)==cund || *(cp+is+1)==cund) icol = 15;
          else {
            cv1 = *cp + (*(cp+1)-*cp)*xx;
            cv2 = *(cp+is) + (*(cp+is+1)-*(cp+is))*xx;
            cv = cv1 + (cv2-cv1)*yy;
            icol = gxshdc(shdlvs,shdcls,shdcnt,cv);
          }
          if (icol!=scol && icol>-1) gxcolr(icol);
          scol = icol;
        }
        uv1 = *up + (*(up+1)-*up)*xx;
        uv2 = *(up+is) + (*(up+is+1)-*(up+is))*xx;
        uv = uv1 + (uv2-uv1)*yy;
        vv1 = *vp + (*(vp+1)-*vp)*xx;
        vv2 = *(vp+is) + (*(vp+is+1)-*(vp+is))*xx;
        vv = vv1 + (vv2-vv1)*yy;
        auv = fabs(uv); avv=fabs(vv);
        if (auv<0.1 && avv<0.1) break;
        if (auv>avv) {
          vv = vv*fact/auv;
          uv = uv*fact/auv;
        } else {
          uv = uv*fact/avv;
          vv = vv*fact/avv;
        }
        x = x + uv;
        y = y + vv;
        ii1 = (int)(x*rscl);
        ij1 = (int)(y*rscl);
        ii1 = ij1*iss + ii1;
        if (ii1<0 || ii1>=siz) break;
        if (*(it+ii1)==1) break;
        if (ii1!=iisav && iisav>-1) *(it+iisav) = 1;
        if (ii1==iisav) iacc++;
        else iacc = 0;
        if (iacc>10) break;
        iisav = ii1;
        // wf ----------------- gxconv_lonlat (x,y,&xx,&yy,SLon, SLat, XInterval , YInterval);//,3);
        if (icol>-1) {
          if (bflg) {gxplot(xold,yold,3); bflg=0;}
          gxplot (xx,yy,2);
        } else bflg = 1;
        xold = xx;
        yold = yy;
        acnt++;
        if (acnt>20) {
          if (icol>-1) strmar (xxsv,yysv,xx,yy);
          acnt = 0;
        }
        xxsv = xx; yysv = yy;
      }
      bflg = 0;
      x = xsav; y = ysav;
      // wf ---------------------------- gxconv_lonlat (x,y,&xx,&yy,SLon, SLat, XInterval , YInterval);//,3);
      gxplot (xx,yy,3);
      xxsv = xx;
      yysv = yy;
      iisav = -999;
      iacc = 0;
      acnt = 19;
      while (x>=0.0 && x<(float)(is-1) && y>=0.0 && y<(float)(js-1)) {
        ii = (int)x;
        jj = (int)y;
        xx = x - (float)ii;
        yy = y - (float)jj;
        up = u + jj*is+ii;
        vp = v + jj*is+ii;
        if (*up==uund || *(up+1)==uund ||
            *(up+is)==uund || *(up+is+1)==uund) break;
        if (*vp==vund || *(vp+1)==vund ||
            *(vp+is)==vund || *(vp+is+1)==vund) break;
        if (flag) {
          cp = c + jj*is+ii;
          if (*cp==cund || *(cp+1)==cund ||
              *(cp+is)==cund || *(cp+is+1)==cund) icol = 15;
          else {
            cv1 = *cp + (*(cp+1)-*cp)*xx;
            cv2 = *(cp+is) + (*(cp+is+1)-*(cp+is))*xx;
            cv = cv1 + (cv2-cv1)*yy;
            icol = gxshdc(shdlvs,shdcls,shdcnt,cv);
          }
          if (icol!=scol && icol>-1) gxcolr(icol);
          scol = icol;
        }
        uv1 = *up + (*(up+1)-*up)*xx;
        uv2 = *(up+is) + (*(up+is+1)-*(up+is))*xx;
        uv = uv1 + (uv2-uv1)*yy;
        vv1 = *vp + (*(vp+1)-*vp)*xx;
        vv2 = *(vp+is) + (*(vp+is+1)-*(vp+is))*xx;
        vv = vv1 + (vv2-vv1)*yy;
        auv = fabs(uv); avv=fabs(vv);
        if (auv<0.1 && avv<0.1) break;
        if (auv>avv) {
          vv = vv*fact/auv;
          uv = uv*fact/auv;
        } else {
          uv = uv*fact/avv;
          vv = vv*fact/avv;
        }
        x = x - uv;
        y = y - vv;
        ii1 = (int)(x*rscl);
        ij1 = (int)(y*rscl);
        ii1 = ij1*iss + ii1;
        if (ii1<0 || ii1>=siz) break;
        if (*(it+ii1)==1) break;
        if (ii1!=iisav && iisav>-1) *(it+iisav) = 1;
        if (ii1==iisav) iacc++;
        else iacc = 0;
        if (iacc>10) break;
        iisav = ii1;
        // wf ------------- gxconv_lonlat (x,y,&xx,&yy,SLon, SLat, XInterval , YInterval);//,3);
        if (icol>-1) {
          if (bflg) {gxplot(xold,yold,3); bflg=0;}
          gxplot (xx,yy,2);
        } else bflg = 1;
        xold = xx;
        yold = yy;
        acnt++;
        if (acnt>20) {
          if (icol>-1) strmar(xx,yy,xxsv,yysv);
          acnt = 0;
        }
        xxsv = xx; yysv = yy;
      }
    }
    i2++;
    if (i2==iss) { i2 = 0; j2++; }
  }
  free (it);

  retArray=write2Array( u, v,  is,  js);

  writeLineOut("teststrm.000f.txt",retArray);
  strmShadFree();
  return retArray;
}
